Shear modulus of the hadron— quark mixed phase 



Nathan K. Johnson-McDanicl^' ^ and Benjamin J. Owcn^ 

^Institute for Gravitation and the Cosmos, Center for Particle and Gravitational Astrophysics, 
Department of Physics, The Pennsylvania State University, University Park, PA 16802, USA 
^ Theoretisch-Physikalisches Institut, Friedrich-Schiller-Universitdt, Max-Wien-Platz 1, 07743 Jena, Germany 

(Dated: July 30, 2012) 

Robust arguments predict tiiat a liadron-quark mixed phase may exist in the cores of some 
"neutron" stars. Such a pliase forms a crystalline lattice with a shear modulus higher than that 
of the crust due to the high density and charge separation, even allowing for the effects of charge 
screening. This may lead to strong continuous gravitational-wave emission from rapidly rotating 
neutron stars and gravitational-wave bursts associated with magnetar flares and pulsar glitches. 
We present the first detailed calculation of the shear modulus of the mixed phase. We describe 
the quark phase using the bag model plus first-order quantum chromodynamics corrections and the 
hadronic phase using relativistic mean-field models with parameters allowed by the most massive 
pulsar. Most of the calculation involves treating the "pasta phases" of the lattice via dimensional 
continuation, and we give a general method for computing dimensionally continued lattice sums 
including the Debye model of charge screening. We compute all the shear components of the elastic 
modulus tensor and angle average them to obtain the effective (scalar) shear modulus for the case 
where the mixed phase is a polycrystal. We include the contributions from changing the cell size, 
which are necessary for the stability of the lower-dimensional portions of the lattice. Stability also 
requires a minimum surface tension, generally tens of MeVfm"^ depending on the equation of state. 
We find that the shear modulus can be a few times 10"^'^ ergcm"^, two orders of magnitude higher 
than the first estimate, over a significant fraction of the maximum mass stable star for certain 
parameter choices. 

PACS numbers: 97.60. Jd, 26.60.Dd, 62.20.de, 04.30.Db 



I. INTRODUCTION 

Asymptotic freedom implies that matter at asymptot- 
ically high densities consists of deconfined quarks, and 
these densities may overlap with the range found in neu- 
tron stars. "'^ This was first noted by Collins and Perry 
[H , although before the discovery of asymptotic freedom 
less concrete suggestions were made by Itoh and Bod- 
mer and the first astrophysically concrete treatment 
(realistic maximum mass, etc.) was given later by Wit- 
ten [1] . More recent developments include the realization 
that any quark matter in neutron stars may be in a color 
superconducting state. The most well-known such state 
is the so-called color-flavor locked (CFL) superconduct- 
ing state present at asymptotically high densities [Bj , but 
there are other color superconducting states that may be 
present at lower densities — see @ for a general review. 
Quark matter could thus have a crystalline structure with 
a very high shear modulus Q , which became relevant to 
gravitational- wave searches several years ago Q . 

Even if none of these mechanisms applies and quark 
matter itself is fluid, there is a generic and robust ar- 
gument that neutron-star cores are filled with a mixed 
phase of hadron and quark matter, which has a crys- 
talline structure due to electrostatic and surface-tension 
effects. Some recent gravitational- wave searches [l^ [T]| 



^ We use the term "neutron star" to refer to any compact star 
made of cold-catalyzed matter, regardless of actual composition. 



have reached sensitivities for which the first rough esti- 
mate of the shear modulus of this mixed phase [l^l — 
which is smaller than the estimates for quark matter — is 
relevant. Therefore it is worth a more careful calculation 
to get a better idea of which gravitational-wave searches 
become relevant to this broader range of theoretical mod- 
els. The justification for the existence of the mixed phase 
is as follows: 

Glendenning [l^ [13 (see also [H ^ for reviews) ar- 
gued that the phase transition from hadrons to quarks 
happens gradually, with a mixed phase existing over a 
wide range of pressures. A neutron star containing such 
a mixed phase is called a hybrid star. The argument for 
this possibility is very robust (provided that the phase 
transition is first order, as is generally expected), relying 
on the fact that charge can be locally separated between 
the two phases while maintaining global neutrality; this is 
energetically favored because it allows the hadron phase 
to reduce its isospin asymmetry without producing lep- 
tons. Observational constraints, even the discovery of a 
1.97 ± 0.04 Mq pulsar do not rule out the possibil- 
ity that large regions of the most massive neutron stars 
are composed of the mixed phase: The input parameters 
of the models have more than enough uncertainty to al- 
low high- mass stars [l8l - l22j . (One can even still obtain 
pure quark cores in such massive stars, with appropriate 
parameter choices [l^, though we do not consider such 
extreme cases here, for simplicity.) 

Like the nuclear "pasta" phases in the crust (first 
discussed by Ravenhall, Pethick, and Wilson [1^), the 
hadron-quark mixed phase is thought to form a lat- 
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tice. This consists of charged blobs of the rare phase 
in a background of the common phase and has a vary- 
ing dimensionahty due to the competition between sur- 
face tension and the Coulomb fo rce, as first suggested by 
Heiselberg, Pethick, and Staubo j2J]. At the lowest pres- 
sure of the phase transition, small quark droplets form 
a 3-dimcnsional lattice in a background of hadrons. At 
higher pressures the droplets grow and give way to a 2- 
dimensional lattice of rods, then a 1-dimensional lattice 
of interleaved quark and hadron slabs. At even higher 
pressures the progression is reversed, with the quark slabs 
outgrowing the hadron slabs, then hadrons forming rods 
and droplets in a quark background before vanishing al- 
together. 

These mixed-phase lattices can have a larger shear 
modulus than that of the crust, due mainly (in three 
dimensions) to the larger charge separation (from several 
hundreds to more than a thousand elementary charges 
per blob, rather than tens for nuclei). The shear mod- 
ulus was first estimated very roughly by one of us fl^ . 
including a simple model of charge screening, the effects 
of which can change the result by orders of magnitude. A 
more detailed calculation was begun by Nayyar [2^ , and 
a rough estimate neglecting charge screening but includ- 
ing the surface energy is given in Sec. 7.7.2 of Haensel, 
Potckhin, and Yakovlcv [26| . The latter also summarizes 
related work on other exotic phases: High shear mod- 
uli for pion condensates have been predicted as early as 
Ref. |27| . and similar estimates were made even earlier 
for solid neutron cores which were proposed to explain 
glitches of the Vela pulsar . 

Due to the high shear modulus of the mixed phase, hy- 
brid stars could sustain much larger deformations than 
normal neutron stars. This has implications not only 
for pulsar glitches but also for gravitational-wave emis- 
sion, both continuous and in bursts associated with mag- 
netar flares [l^. Upper limits on the gravitational- 
wave energy emitted in magnetar flares have entered the 
range of theoretical predictions (up to 10^^ erg for hybrid 
stars) [n], m, |30|- And the most interesting upper 
limits on continuous-wave emission, those that beat the 
indirect limits (see [Sll, [s^] for reviews), are for several 
stars [1, within an order of magnitude of the 10~^ 

maximum ellipticity first estimated for hybrid stars (l2j 
and at the level of the 10^^ maximum that applies if more 
recent results on the maximum crust breaking strain (33 | 
hold for the mixed phase. Therefore, more detailed calcu- 
lations of the crystalline structure of the mixed phase are 
interesting for the interpretation of gravitational-wave 
observations even now, and will become more so when 
the upgraded "advanced" detectors come on-line. 

Here we improve upon previous estimates of the mixed- 
phase shear modulus O, IH, [1^ with the first detailed 
calculation. The implications for magnetar flares were 
discussed by Corsi and Owen (ssf . We will present an 
improved calculation of the consequences for continuous 
waves elsewhere (36j . 

We use the standard models used by Glendcnning [l^ 



for the hadronic and quark equations of state (EOSs), 
as well as the standard Gibbs method for calculating 
the bulk properties and lattice structure of the mixed 
phase. We thus use a relativistic mean field model for 
the hadronic matter. As Norsen and Reddy [s^l note, 
this assumption can be of dubious validity for small re- 
gions of hadronic matter, since the fields do not take on 
their mean values in that case. However, this is not much 
of a concern for us, since, as Norsen and Reddy discuss, 
it is only a large effect for lower-dimensional hadronic 
blobs, and these are only present in stable stars for one 
of the EOSs we consider, where they have relatively large 
dimensions. (See Sec. |TT] for further discussion.) We do 
not have to worry about such errors in our description 
of the quark blobs, since we are using an improved bag 
model for them. We also include the oft-neglected contri- 
bution from the surface tension to the pressure balance 
in a few of our EOSs, following [STi - lioj . 

We obtain significantly larger shear moduli than 
Ref. [l2j primarily due to our inclusion of the contri- 
butions that arise from the change in the cell volume 
when one shears the lower-dimensional lattices. (Some 
of these contributions were considered by Pethick and 
Potekhin [4l| for the case of the nuclear pasta in the 
crust. Additionally, the estimate of Haensel, Potekhin, 
and Yakovlev [2y] is based on these contributions.) 
This is the most important improvement in our anal- 
ysis. These contributions are necessary for the lower- 
dimensional lattices to be stable, and this stability leads 
to an EOS-parameter-dcpendcnt minimum surface ten- 
sion. Additionally, for large but reasonable surface ten- 
sions, these contributions significantly increase the over- 
all shear modulus of the lower-dimensional lattices. 

Another fundamental improvement is that we treat the 
varying dimensionality of the pasta phases using dimen- 
sional continuation. This technique was first suggested 
by Ravenhall, Pethick, and Wilson [2^ for the nuclear 
pasta in the crust, and then applied to the hadron-quark 
pasta by Glendcnning (l5| . It is a relatively simple way 
to approximate the complicated structure seen in, e.g ., 
molecular dynamics simulations of nuclear pasta |42l445l | . 
In order to perform the dimensionally continued lattice 
sums, we have developed a generalization of the standard 
Ewald method [1^ (see [43] for a modern treatment). 
The Ewald method was used by Fuchs [i^ in his pio- 
neering calculation of clastic coefficients, and has very 
recently been used by Baiko [l^ in a calculation of the 
effective shear modulus of the neutron star crust. The 
underlying dimensionally continued Poisson summation 
formula has been reported elsewhere by one of us jsof . 
Here we describe the practical implementation to the 
situation at hand, including the details of the "Ewald 
screening function" and the dimensional continuation of 
the specific family of lattices we consider. 

We would have liked to dimensionally continue the 
family of root lattices, AJ, which solve the covering prob- 
lem in dimensions d < 5 (see, e.g.. Sec. 1.5 in Conway and 
Sloane [EH). However, as discussed in Sec. IV Bl there is 
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no obvious way of doing so. Instead, we split the lattice 
up into liyperlattices and introduce a freely specifiable 
interpolation function for the hyperlattices' separation. 
We find that our results are quite insensitive to the spe- 
cific choice of this function (particularly given the much 
larger uncertainties in the input parameters). 

In yet another improvement, we perform rather than 
guess the angle averaging to obtain a scalar shear mod- 
ulus from the (anisotropic) elastic modulus tensor. This 
makes the assumption (standard in the literature) that 
there are many regions with random orientations, so that 
the mixed phase can be treated as a polycrystal. The 
magnetic fields present in a neutron star may cause this 
assumption to be violated, and relaxing it would be an 
interesting topic for further investigation. 

We also use the full Debye (linear) model for charge 
screening rather than multiply the unscreened result by 
a simple correction factor. Even this is a simplified treat- 
ment of screening effects, particularly because it is lin- 
ear. Also, because we treat the blobs as point charges, it 
means that we do not include the contribution of charge 
screening in our computation of the blobs' energy (which 
is used to obtain the blob size and lattice spacing). How- 
ever a more detailed treatment (following, e.g., j37H4C| ) 
would be much more computationally intensive than the 
remainder of our approach, so we have left such inves- 
tigations to future work. (Recent work by Endo (ioj 
shows that the mixed phase can still occupy much of the 
star even with a nonlinear treatment of charge screening, 
though one cannot draw any definite conclusions about 
our results from Endo's because he uses different EOS 
parameters.) Some indication of the magnitude of the 
error we make in using the Debye model can be seen in 
the jumps in the shear modulus in the figures in Sec. I VII 

We use electrostatic units and set /i = c = 1, so we 
will generally express masses in MeV and lengths in fm. 
All the computations were performed using Mathemat- 
ICA 7's default methods to solve algebraic equations, nu- 
merically evaluate integrals, etc. 

The paper is structured as follows: In Sec.|ll]we review 
the models we use for the hadronic and quark EOSs, 
as well as how we compute the bulk properties, lattice 
structure, and charge screening of the mixed phase. In 
Sec. lIIIi we give an overview of the elastic response of the 
lattice for the various integer dimensions, and show how 
we compute an average shear modulus in a dimension- 
ally continued manner. We describe how to compute the 
relevant elastic constants in Sec. IIVI and how to dimen- 
sionally continue the resulting lattice sums in Sec. [V] We 
give our results for the shear modulus in Sec. IVIi along 
with some discussion, and conclude in Sec. I VIII We dis- 
cuss various checks on our computations of the lattice 
sums in the Appendix. 



II. MIXED PHASE 

We model the hadronic and quark phases following 
Glendenning [3, using a relativistic mean field the- 
ory model for the hadronic matter (from Glendenning's 
Chap. 4), and an improved bag model description of the 
quark matter (from Glendenning's Chap. 8). (For sim- 
plicity, we do not consider the possibility of color su- 
perconductivity.) We have chosen to compute the EOSs 
"from scratch" instead of using any sort of tabulated EOS 
(except for the low-density EOS, which has a negligible 
effect on our results). This allows us to include the effects 
of surface tension on the EOS (as discussed in Sec. Ill C|) . 
and to investigate the effects of different EOS parame- 
ters on the shear modulus. It also lets us compute the 
Debye screening length, for which we need to know how 
the chemical potentials of the particle species vary with 
density. The models we use for the hadronic and quark 
phases are relatively simple, but should contain at least a 
rough description of the relevant physics for the situation 
under consideration. And given the significant uncertain- 
ties associated with any description of cold, dense matter, 
even relatively simple models should allow an adequate 
sampling of the relevant parameter space. 

As mentioned previously — and emphasized by Norsen 
and Reddy [s^] — the mean field theory description is not 
accurate for small regions of hadronic matter. However, 
this is of little concern for us, since this should only be a 
large effect for small, low-dimensional hadronic blobs, as 
discussed by Norsen and Reddy. The only EOS for which 
one obtains such blobs in stable stars, LKRl, has low- 
dimensional hadronic blob radii of ^ 10 fm, as illustrated 
in Sec. Ill CI Taking the hadron-kaon results from Norsen 
and Reddy to apply to our case, we see that the largest 
finite-size corrections will thus be ^ 10%. 

These EOSs depend on a variety of parameters, dis- 
cussed in more detail in the following subsections. We 
consider a small collection of representative sets of pa- 
rameters given in Table HI All of these parameter sets are 
chosen to yield a maximum Oppenheimer-Volkov (OV) 
mass compatible with the recent Demorest et al. [l3l mea- 
surement of a 1.97 ± O.O4M0 neutron star. 

Since the hadronic EOS we consider is only valid at 
densities well above neutron drip, we add on a stan- 
dard low-density EOS for baryon number densities ub < 
0.08 fm~'^. This is the combination of the Baym, Pethick, 
and Sutherland (BPS) ^ EOS for ub < 0.001 fm"^ 
and the Negele and Vautherin $^ EOS for 0.001 fm"^ < 
Ub < 0.08 fm~'^ used by Lattimer and Prakash 1541 . We 
use the table provided by Kurkela et al. [s^ at j56l |. The 
precise choice of the low-density part of the EOS has a 
negligible effect on the maximum mass of a stable star. 

Since we are interested in how large astrophysical ef- 
fects could be, we generally choose our EOS parameters 
to be middle-ground estimates, but also consider some 
more extreme cases, to yield massive stars with a large 
region of mixed phase. (However, we do not consider 
stars with quark cores, for simplicity.) The Hyl EOS is 



4 



taken from Nayyar [25| . correcting an error in his code 
which led to a lower maximum mass. The Hyl' EOS 
changes the quark bag constant to obtain a larger region 
of mixed phase with more pasta phases in stable stars. 
The Hyl/i and Hylcr EOSs each change the treatment 
of one portion of the calculation to give some indica- 
tion of how much these affect our results: Hyl/x uses a 
chcmical-potential-depcndent quantum chromodynamics 
(QCD) renormalization scale, while Hylcr includes the 
effects of surface tension on the pressure balance of the 
two phases; Hyl/ia includes both. 

Following Weissenborn et al. [23|, we also consider a 
case called LKRl with the NL3 hadronic EOS parameters 
of Lalazissis, Konig, and Ring (LKR) [s^]. These param- 
eters give a very stiff hadronic EOS (a purely hadronic 
star would have a maximum mass of 2.78M0) and thus 
allow for neutron stars with a large region of mixed phase 
that still are compatible with a 1.97 ± 0.04 neutron 
star. For the corresponding quark EOS we picked pa- 
rameters that lead to stable stars that include all the 
pasta phases. We also considered a case with somewhat 
generic parameters, generally picking them to be around 
the midpoint of the accepted range — this is the "generic" 
parameter set. These were not fine-tuned at all. Ad- 
ditionally, we have considered a variant that uses these 
more modern hadronic EOS parameters along with a bag 
constant and QCD coupling constant that yield a larger 
region of mixed phase — this is the generic' parameter set. 

We plot these EOSs up to the maximum density 
present in a stable OV star in Fig. [T] We do not show the 
various flavors of the Hyl EOS here, since the resulting 
traces are all very similar, though we do show the small 
violations of le Chatelier's principle for the Hylcr and 
Hyl^CT EOSs, compared with the Hyl EOS in Fig. H 
(We do not show the trace for the Hyl/x EOS, as it is 
almost identical to that for the Hyl EOS.) 
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FIG. 1: Pressure vs. energy density for the EOSs from Tabled 
plotted up to the maximum energy density present in a stable 
OV star. 
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FIG. 2: Pressure vs. energy density for the Hyl EOS and two 
flavors that include the surface tension contribution to the 
pressure balance, illustrating the violations of le Chatelier's 
principle. (Again, each of these is plotted up to the maximum 
energy density present in a stable OV star.) 



A. Hadronic EOS 

We construct the hadronic EOS following the recipe 
described in Chap. 4 of Glendenning's book [la] through 
Sec. 4.9. We thus use a relativistic mean field description, 
with a Lagrangian that contains the standard scalar and 
vector (cr and u)) fields, as well as scalar self-interactions, 
and the isospin asymmetry force, mediated by the vec- 
tor meson p. We include neutrons, protons, electrons, 
and muons. See [2^ and [11] for further details about 
the general framework and calculational procedure we 
used. However, we have used updated parameters in the 
models, as discussed above. Note that these references 
also consider hyperons, which we do not include in the 
EOSs considered here. In particular, hyperons and the 
hadron-quark mixed phase have been found to be mutu- 
ally exclusive in some studies (see, e.g., pTI. [59j). 

The model has five input parameters: Two are known 
reasonably well, viz., the number density and binding en- 
ergy per nucleon of nuclear matter at saturation, no = 
0.16 ± 0.01 fm"^ and {EjA)^ = -16 ± 1 MeV. The 
remaining three are not nearly so well known: The nu- 
clear incompressibility K is thought to lie between 200 
and 300 MeV, with many authors placing it in the range 
240 ±10 MeV (see [13]), while the (scaled) Dirac effective 
mass of the nucleons. m* /m, is thought to be between 
- 0.53 and 0.96. (Here m = 938.93 MeV = 4.7582 fm"^ 
is the average of the neutron and proton masses.) The 
symmetry energy Osyrn is thought to lie between 28 and 
34 MeV (see Li et aL |6l|). 

All parameter ranges are taken from Eq. (88) in 
Steiner et al. [i^ unless otherwise noted. The range for 
the Dirac effective mass is computed from that for the 
Landau effective mass given by Steiner et al., following 
Glendenning's Eq. (4.117) [l6j . with a Fermi momentum 
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TABLE I: EOS parameters and properties of their associated OV stars. In the Aqcd column, fig denotes the cases where 
we take the QCD renormalization scale to be given by the average quark chemical potential at each density, as discussed in 
Sec. Ill Bl In the a column, we have denoted the cases in which we include the surface tension contribution to the pressure 
balance by a "p." M^'j,-^^'"* gives the masses of stars that first contain hybrid matter; -R|5ax"*/^ denotes the maximum radius 
fraction occupied by hybrid matter (i.e., the radius fraction for the maximum mass star); and Cmax denotes the maximum 
compactness (2GM / R(?) of a star. We also give the composition of the rare phase ("Q" for quark and "H" for hadronic) and 
the dimension of the lattice at the center of the maximum mass star. See the text for the definitions of other parameters. 



at saturation of 1.33 ± 0.03 fm""'^, from the Steiner et al. 
result and Glendenning's Eq. (4.110). 

Note that almost all of the NL3 parameters from LKR 
fall outside of the ranges we have given. We still con- 
sider them for two reasons. First, the LKR NL3 EOS is 
often treated as the paradigmatic stiff hadronic EOS in 
the literature; in particular, we were inspired to consider 
these parameters by the LKR EOS's recent use in Weis- 
senborn et al. Second, the LKR NL3 parameters are 
still regarded as providing an excellent fit to the nuclear 
binding energy and charge radius, explaining why this 
EOS remains in use. 



B. Quark EOS 

The quark matter is described by the MIT bag model 
(first given in Chodos et al. [11], and discussed briefly 
in Glcndcnning [T6j). to capture the basic physics of con- 
finement. This is supplemented with first-order QCD cor- 
rections to the thermodynamic potential for free quarks 
from Far hi and JafFe [6J], as described in Chap. 8 of Glen- 
denning [iBl-^ We take all the quarks to be massive (as 
opposed to the usual treatment in which onl y th e strange 
quark is taken to be massive — see, e.g., |20|. l65l [66j — but 
note that Christiansen and Glendenning |67| also take 
all three quarks to be massive). Sec Chap. 6 in [2^ for 
further details of the calculation. 

Here the parameters are all quite uncertain. The only 
firm constraint is that nonstrange quark matter is not 
absolutely stable, since this contradicts the observed ex- 
istence of nuclei composed of nucleons. In addition, we 
assume that strange quark matter is not absolutely sta- 
ble, since we consider hybrid stars, not strange stars. In 



the pure bag model (with no QCD corrections — i.e., with 
zero QCD coupling constant), this implies that the fourth 
root of the bag constant, B'^/^, is greater than 145 MeV. 
When one includes QCD corrections, the minimum bag 
constant decreases — sec, e.g., the discussion in Alford et 
al. |65j . There is no known upper bound on though 
for sufficiently large values of B (with all other EOS pa- 
rameters fixed), stable stars do not contain deconfined 
quark matter. 

We also need to know the QCD coupling constant and 
quark masses at the energy scale Aqcd at which we renor- 
malize the perturbative contributions to the thermody- 
namic potential. This energy scale is typically taken 
to be the scale of the quark chemical potentials in the 
problem. Glendenning uses 300 MeV, though we find 
that the chemical potentials in the situations we con- 
sider are somewhat higher (up to ^ 480 MeV in the 
densest regions of stable stars). Therefore we also con- 
sider a density-dependent renormalization scale, given by 
Aqcd = ■= iti-p + fJ'n)/6, where fXp and /i„ are the 
proton and neutron chemical potentials, respectively, so 
flq is an average quark chemical potential. (This was 
inspired by the chemical potential-dependent renormal- 
ization scale used in [6^, though the specifics of their 
treatment differ.) 

At this scale, the QCD coupling constant as is large 
and not well known (see Fig. 1 in [66| for the results from 
running the coupling constant using relatively low-order 
beta function expressions). We have considered values 
between 0.5 and 0.7, which are quite middle-of-the-road 
(particularly compared with the naive beta function re- 
sult of much greater than 1, though there is evidence 
that the QCD coupling constant "freezes" to a value of 
less than 1 at low energies — see, e.g., [6^ H^)- Weis- 
senborn et al. [201 consider values from to 0.94.'^ while 



^ Note that Glendenning corrects a sign error in Farhi and JafTe's 
result for the thermodynamic potential in his Eq. (8.14). 



^ This is converted from Wcisscnborn et al.'s range of 0.4 to 1 
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Nayyar [2^ uses values of 0.45 and 0.6. The quark masses 
are similarly uncertain. There are significant uncertain- 
ties (up to ~ 50%) even at the 2 GeV scale at which 
the Particle Data Group [t^ quotes the masses, and ad- 
ditional uncertainties in running them to lower energies. 
For simplicity, we have taken the up, down, and strange 
quarks to have masses of 2.5, 5, and 100 MeV, respec- 
tively, around the median of ranges for the MS current 
masses at 2 GeV given in the 2010 Review of Particle 
Physics [zOl • We also consider a few cases with a strange 
quark mass of 150 MeV, for continuity with Nayyar [25| . 



C. Hybrid phase and its lattice structure 

We determine the bulk properties of the mixed phase 
using the Gibbs equilibrium conditions, following Glen- 
denning (see, e.g., Chap. 9 in [3): The appearance of 
the mixed phase at a certain baryon number density is 
signaled by pure quark matter having a larger pressure 
than hadronic matter. We then determine the volume 
fractions of the two coexisting phases by demanding the 
equality of the phases' pressures and chemical potentials, 
in addition to global electrical neutrality. We compute 
the mixed-phase bulk pressure, energy, and baryon densi- 
ties by weighting the contributions from each phase using 
their volume fractions. 

We also need to obtain the crystalline structure of 
the mixed phase. This requires knowing the surface 
tension a of the hadron-quark interface. While a has 
been estimated very roughly in the literature (e.g., [13] 
and the references given in Endo its value is still 

quite uncertain, so we simply take it to range from 10 
to 80 MeV fm~^. (The lower bound is the default value 
used by Nayyar [25|, and the upper bound is the default 
value used by Glendenning We make sure that it 

does not make an appreciable contribution to the overall 
energy density, since, as Glendenning discusses, the sys- 
tem's energy should not increase upon opening up a new 
degree of freedom. In the cases we have considered, the 
blobs' energy density is < 2% of the total. Similarly, the 
lattice's electrostatic pressure is also < 3% of the bulk 
pressure. 

Although the Coloumb and surface energy of the blobs 
is only a small contribution to the overall energy density, 
it can still be significant compared to the energy differ- 
ence between the pure hadronic phase and the mixed 
phase at a given baryon density, particularly at lower 
densities. Thus, as discussed by Heiselberg, Pethick, and 
Staubo [H] and Alford et al. [7li, for sufficiently high 
surface tensions, the mixed phase may be disfavored com- 
pared to the sharp transition predicted by the Maxwell 



for Alford et ai.'s [65^ using the massless quark expression of 
Os = (1 - 04)77/2 [cf. Eq. (3) in Alford et al. [H and Eq. (6.4) 
in Nayyar 1251 ]. 



construction. Christiansen and Glendenning, however, 
argue that the Maxwell construction should always cor- 
respond to an excited state for a multi-component system 
such as we are considering, and if any models predict oth- 
erwise for certain parameter values, they are not accurate 
for those parameters [63, [zS] ■ 

Regardless, these energy arguments are all local (i.e., 
at a fixed baryon density) . If one considers the energy of 
the star as a whole (at a fixed total baryon number) , then 
the mixed phase can still be favored, particularly if one 
accounts for the increase of the star's own gravitational 
binding energy due to the attendant softening of the EOS 
at high densities. We shall present our analysis of the ex- 
tent to which the mixed phase is present in stable stars 
for various surface tensions when we consider the associ- 
ated maximum quadrupoles in [36j . Suffice it to say that 
given all the uncertainties present in these calculations, 
we still feel comfortable quoting results for surface ten- 
sions of 80 MeV fm~^, at least as an indication of the 
upper bounds possible for the shear modulus. (We also 
discuss the dependence of the shear modulus on surface 
tension in Sec. IVII ) Note that if one includes the lattice 
contributions to the EOS (energy density and pressure, 
plus the energy density of the blobs themselves) , one ob- 
tains surface tension-dependent corrections to the values 
for stellar properties given in Table H] However, these are 
at most a few percent for the surface tensions we consider 
here. 

Returning to the lattice structure, we note that the 
competition between the aforementioned surface tension 
and the Coulomb energy of the charges, the charge breaks 
up into blobs, which will then form a lattice. In order to 
determine the properties of this lattice, we employ the 
Wigner-Seitz approximation, which divides the lattice 
up into noninteracting, electrically neutral cells (approxi- 
mating the lattice's Voronoi cell by a sphere of equivalent 
volume). Each cell contains a charged blob at its center 
surrounded by an equal amount of compensating charge. 
At lower densities, the blob is formed out of quark mat- 
ter, and the hadronic matter provides the compensating 
charge; at higher densities — not reached in stable stars 
for many EOSs — the roles are reversed. The ratios of 
the volumes of the blob and the cell are fixed, since the 
baryon and quark volume fractions are fixed. If we de- 
note the volume fraction of the rare phase by x, then we 
have X = {r/RY, where r and R are the radii of the blobs 
and cells, respectively (the half-thickness of the slabs in 
the one-dimensional case), and d is the dimensionality of 
the cell and lattice. If we denote the volume fraction of 
the quark phase by Xi we have 

\l~X, X> 1/2. ^ ' 

We then determine r and d by minimizing the cell's en- 
ergy per unit volume. 

Following the suggestion initially made by Ravenhall, 
Pethick, and Wilson [2^ in the nuclear pasta case, we 
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take the lattice's dimension to be a continuous variable. 
In the integer dimension cases, we have A^: a body- 
centered cubic (bcc) lattice of spherical drops in 3 di- 
mensions, a hexagonal lattice of cylindrical rods in 2 di- 
mensions, and equally spaced rectangular slabs in 1 di- 
mension. (These lattices are seen in the molecular dy- 
namics simulations of nuclear pasta [4^ - |45j . In partic- 
ular, Watanabc et al. (4^] obtain a hexagonal lattice of 
rods by adiabatically compressing a bcc lattice of drops 
in their molecular dynamics simulations.) 

The expression for the energy per unit volume of a 
cell is then given by [Eqs. (9.19), (9.23), and (9.24) in 
Glendenning [l^ ] 



(2) 



where Q is the cell volume, and C{x) and S{x) correspond 
to the Coulomb and surface contributions, respectively, 
with 

Cix) := 2^[(te - qQ){x)?xh{x), S{x) := xad. (3) 

Here qn and qq arc the hadron and quark charge densi- 
ties and 



fd{x) :-- 



1 



d+2 



2 - dx^-^'^ 
d~2 



(4) 



The singularity at d — > 2 is removable, since the limit ex- 
ists [and has the expected value; sec Eq. (9.30) in Glen- 
denning [l6j ]. One can now immediately read off the 
minimizing r [Eq. (9.27) in Glendenning jla | ], viz.. 



Six) 



-,1/3 



2C{x) 



(5) 



One then obtains d S [1,3] by minimizing (over 
that range) the resulting expression for the cell en- 
ergy (using the minimizing r), viz., Ecdi/^ = 
(3/2)[2C'(x)]i/3[S'(a;)]2/3 [see Eq. (9.28) in Glenden- 
ning [16| for an explicit expression in terms of x, etc.]. We 
show the dependence of blob radius and lattice spacing 
[given in Eq. (|49p ] on x for a few representative EOSs in 
Fig. 131 (The blob radii and lattice spacings of the three 
EOSs shown span the range covered by the entire set of 
EOSs in Table HI) 

As discussed in, e.g., [37l - l40j . the surface tension con- 
tributes to the pressure equality. Explicitly, we have a 
pressure difference between the dominant and rare phases 
of [cf. Eq. (A5) in ^ ] 



^dominant Pr. 



(d- 1)ct 



(6) 



where r is the radius of the blobs. We have included this 
in our treatment of certain EOSs, though in a somewhat 
simplified version, so that it is numerically tractable. One 
obtains these EOSs by solving the appropriate equations 
at progressively larger and larger baryon densities (as 
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FIG. 3: The blob radius and lattice spacing versus x for three 
representative EOSs. 



discussed by Glendenning |l6|). Our simplification con- 
sists of taking the lattice's dimension at a given baryon 
density to be fixed at the value obtained at the previ- 
ous baryon density (instead of solving for the dimension 
including the dimension-dependent surface tension con- 
tribution to the pressure equality). Since this correction 
to the pressure equality is largest at the lowest-density 
portions of the mixed phase (i.e., the quark drop por- 
tions, with d = 3), our procedure for d seems likely to 
account for the primary effects from the surface tension's 
contribution to the pressure equality. (Note that we take 
the bulk pressure to be that of the dominant phase.) 

With this treatment of the surface tension contribution 
to the pressure equality, we find a slight violation of le 
Chatelier's principle (i.e., monotonic increase of pressure 
with energy density) close to the hybrid transition. This 
is not present if one does not include the surface ten- 
sion contribution to the pressure equality (see Fig. [5]). 

nations — 
23 — does 



Note that our method of solving the OV ec 
using the enthalpy form given by Lindblom 
not allow for such violations of le Chatelier's principle, 
since it requires one to express the energy density as a 
function of the pressure. We nonetheless quote results 
using the "smoothed-out" version of the EOS produced 
by Mathematica's interpolation, feeling that they still 
give a reasonably accurate depiction of the star's bulk 
properties, since the violations of le Chatelier's principle 
we are considering are small. 



D. Charge screening 

The previous discussion has assumed that the charge 
is uniformly distributed within each blob and in the neu- 
tralizing background. This is not the case, in practice: 
The minimum-energy configuration will have a nonuni- 
form charge distribution, as is discussed in, e.g., (37l - [39j . 
This nonuniform charge distribution — often treated in 
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the perturbative regime as charge screening — will affect 
the cell energy (and thus the lattice properties, for a given 
energy density), as well as the lattice's electrostatic en- 
ergy, and both of these affect the shear modulus. Here we 
only consider the effects on the lattice's electrostatic en- 
ergy due to linear charge screening of point charges, using 
a screened potential. As is commonly done in treatments 
of the mixed phase (e.g., those by Glcndcnning) , we do 
not concern ourselves at all with the rearrangement of 
charge inside the blobs, though this is likely a significant 
effect (as is discussed in |37l - [39| ). 

In this linearized version, we treat screening as a small 
perturbation on the overall energy, leading to the expres- 
sion for the Debye length given in Eq. (1) of [l^l, viz.. 



A = 



dUa 



1 -1/2 



(7) 



Here Qa, n^, and are the charge, number density, and 
chemical potential of particle species a, and the partial 
derivative is evaluated holding the chemical potentials of 
all species besides the ath fixed. 

To include the effects of charge screening on the elec- 
trostatic energy, we dimensionally continue the standard 
screened potential equation (Yukawa in three dimen- 
sions), viz.. 



(Ad-A-2)0 = _4^g5W, 



(8) 



where is the d-dimensional Laplacian, Q is the charge 
of the blob, and 5^'^^ is the d-dimensional Dirac delta 
distribution. 

The resulting potential is 



(r) 



2Q 



(27rAr)'^/2 



-Kd/2- 



(9) 



where is the modified Bessel function of the second 
kind of order v. This expression can be obtained most 
easily by noting that it is the same as the Euclidean scalar 
propagator (up to overall factors), which is g iven in a 
dimensionally continued form in Eq. (5) of |74{.^ 

This potential clearly represents a simplified treat- 
ment of charge screening even in the perturbative limit, 
since we have used a point charge, instead of the ex- 
tended charge distribution of a realistic blob with differ- 
ent screening lengths inside and outside. (We similarly 
use point charges when calculating the potential energy 
of other blobs due to this screened potential.) However, 
it is a practical way for us to account for some of the 
effects of charge screening. Since we are using this point 
charge approximation, we compute the screening length 
using the leptons and the background phase. This is be- 
cause the rearrangement of charge within the blob has no 



* When comparing the two expressions, note that K- 



effect on the blob's electrostatic potential except through 
the effects of screening on the blob's size, which we are 
ignoring here. This switch at x = 1/2 causes a significant 
jump in the screening length there. This translates into a 
jump in the elastic constants, and thus the effective shear 
modulus, as seen in the figures in Sec. lVIl The magnitude 
of this jump gives some indication of the overall error we 
incur through our simplified treatment of charge screen- 
ing in the potential. (However, this likely does not give 
any indication of the errors incurred by neglecting the 
contribution of charge screening to the cell energy, which 
we shall see affects the effective shear modulus both di- 
rectly and indirectly.) 



III. ELASTICITY 

The elastic response of a crystalline lattice is described 
by its elastic modulus tensor. Skips, defined by 

£ = fo + SkjUkl + ^SklpsUklUps, (10) 

where £ is the deformed lattice's energy density, Eq is 
its undeformed energy density, Uki :— dkSxi is the dis- 
placement gradient (where Sxi denotes the displacement 
field), and Ski is the first-order piece of the expansion. 
We have a nonzero first-order piece here, since the un- 
deformed lattice has a nonzero pressure — see Eq. (7) in 
Baiko [4§|. (Wallace [zl] gives further discussion.) Ex- 
plicitly, we have Ski = —Pcs^ki, since the pressure is 
isotropic. (Note that we shall only consider the contri- 
bution from the lattice's electrostatic pressure Pes hi our 
discussion.) Now, we are only interested in the shear 
stress the lattice generates in response to a shear defor- 
mation, and this is given by a different elastic modulus 
tensor, viz., 



Bkips = Skips — Pcs{SksSip — SkiSps). 



(11) 



Specifically, the stress generated by a deformation Uki 
is given by Bki„s(u„s + Usp)/2. [See Eqs. (29) and (30) 
in Baiko [4^.] We can write the components of Bkips 
more simply by using Voigt notation to map them to the 
two-index object Ca/3 using the index mapping given by 
{xx,yy,zz,xy,xz,yz} O {1,2,3,4,5,6} where {x,y,z) 
are Cartesian coordinates. With this notation, we have 
cii = S'liii and C44 = 51212, while C12 = 5*11; 



= K„ 



[Note that Baiko's Cki denotes a different tensor.] 

Since we are interested solely in the lattice's response 
to shears, we can focus our attention on a few compo- 
nents of this tensor. If we just consider the cases where 
a shear strain yields a proportional shear stress (as is the 
case for isotropic materials), then we will have contribu- 
tions from the simple shear portions of the tensor, viz., 
C44, C55, and C66, in addition to the elongational shears 
A12 := (cii + C22)/2 - C12, ^13 := (cn -|- 033)72 - C13, 
and A23 := (C22 + C33)/2 — C23. The simple shears corre- 
spond to stresses or strains of the form 2x^kyi) where Xk 
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and yk are unit vectors in the a;- and y-directions, and 
parentheses on the indices indicate symmetrization. The 
elongational shears correspond to stresses or strains of 
the form XkXi — ykVi — i-e., rotations of the simple shears 
by 7r/4. In the past, some investigations of shear mod- 
ulus effects in the crusts of neutron stars have used the 
simple shear portion of the tensor, C44, as the shear mod- 
ulus of a bcc lattice (for which C44 = C55 = cge) — see 
the discussion in Sec. 7.1 of Chamel and Haensel [tg} . 
As Strohmayer et al. [t^] emphasize, it is inappropriate 
to simply use one component of the elastic modulus ten- 
sor for this purpose. While a truly detailed calculation 
would use the full elastic modulus tensor, if one consid- 
ers a polycrystal consisting of many randomly oriented 
domains, it is possible to use an angle-averaged version 
of the shear portions of this tensor and obtain an upper 
bound. 

This upper bound is due to Voigt [zl], and involves 
the elastic constants given above — see Hill [t^ for the 
proof that the Voigt expression gives an upper bound. 
We do not consider any of the more involved sharper 
bounds, such as the oft-used ones by Hashin and Shtrik- 
man [sO]. (See [8l| for a review of such bounds, and [s^ 
for a Hashin- Shtrikman bound for orthorhombic crystals, 
such as the ones we consider here.) The Voigt result was 
rediscovered by Ogata and Ichimaru [83| in considering 
the shear modulus of the neutron star crust. (They used 
it as an average, as was originally proposed by Voigt, 
not as an upper bound, as in Hill.) Following Ogata and 
Ichimaru [83| and most subsequent work, we shall use the 
Voigt average for our effective shear modulus, giving 

Aicff = {A12 + + yl23)/15 + (C44 + C55 + c66)/5. (12) 

We note, as docs Baiko [i^, that the lattice will likely 
tend to align itself with the star's magnetic field, leading 
to large-scale ordered structure. However the standard 
averaging approach we use here is a necessary first step, 
and may be a good approximation to bulk properties of 
the star if the internal magnetic field lines are tangled as 
in some simulations such as Braithwaite and Spruit [s^ . 

To dimcnsionally continue this expression, we need to 
consider the effects of different lattice dimensionalities on 
the shear elastic coefficients. We shall first discuss the 
integer dimension cases, and then present the dimension- 
ally continued expression for the effective shear modulus 
at the end. (See Pethick and Potekhin 4l[ for related 
discussion about the elasticity of liquid crystals as ap- 
plied to the pasta phases in the crust.) Our discussion 
will be simplified by the fact that we can take the elastic 
constants of all the lattices to have cubic symmetry — 
this holds for the integer dimension lattices and remains 
true for the dimcnsionally continued ones by fiat. How- 
ever, there will be differences between the elongational 
shear elastic constants in which the shear takes place 
completely within the lattice and those that involve an 
elongation along the rods or slabs. We also will find that 
the number of nonzero shear elastic constants is reduced 
due to the translational symmetry along the rods or slabs 



in lower dimensions. Specifically, 

1. In the three-dimensional case, all the shear elastic 
constants are nonzero, and, by cubic symmetry, the 
simple and elongational shear elastic constants are 
each all equal (i.e., C44 = C55 = cge and A12 = 

^13 = ^23 =: Aat)- 

2. In the two-dimensional case, only one of the sim- 
ple shear elastic constants is nonzero, since simple 
shears along the rods do not lead to a shear stress — 
i.e., we have cgg 7^ C44 = C55 = 0, where we have 
taken the rods to point in the z-direction. The 
elongational shears are all nonzero, but are not all 
equal: In general, the one elongational shear per- 
pendicular to the rods (i.e., within the lattice, so 
denoted ^lat) has a different elastic constant than 
do the two elongational shears that involve elon- 
gations along the rods (i.e., perpendicular to the 
lattice, and so denoted A±_). Explicitly, we have 

Alat A12 ^ ^13 = ^23 =: A^. 

3. In the one-dimensional case, all of the simple shear 
elastic constants are zero, due to the translational 
symmetry along the slabs, and the only nonzero 
elongational shear constants are those that involve 
elongation along the slabs. (We neglect any change 
in energy due to shearing one of the slabs.) We 
thus have C44 = C55 = cgg = and, taking the x- 
direction to be perpendicular to the slabs, ^lat := 
A12 = Aiz ^ A23 = 0. 

We can now obtain the dimcnsionally continued ver- 
sion by using the dimcnsionally continued versions of 
the relevant basic combinatorial results: The number of 
independent shears (either simple or elongational) com- 
pletely within the lattice (i.e., perpendicular to the rods 
and slabs in the integer dimension cases) is given by 
d{d— l)/2; this gives the multiplicity of the contributions 
to /ioff- from C44 and Ajat. The number of independent 
elongational shears with one elongation perpendicular to 
the lattice and one within it is given by d(3 — d)\ this 
gives the multiplicity of the contributions from A±_ . The 
dimcnsionally continued version of Eq. (|12p is thus 



d 
15 



1 



Aut + (3 - d)A, 



d{d-l) 
10 



C44. (13) 



IV. CALCULATION OF THE ELASTIC 
CONSTANTS 

Now looking at how to compute Ajat, and C44, 
we can follow Fuchs [48| in expressing these quantities in 
terms of lattice sums. Specifically, we write the electro- 
static energy of an appropriately deformed lattice as a 
lattice sum and takes derivatives with respect to the de- 
formation to obtain lattice sum expressions for the elastic 
constants. In Ax, we will also need to consider the con- 
tributions from the cell energy [obtained using Eq. 
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One also expects there to be contributions to Ai at an d C44 
due to changing the shape of the unit ceh [63, 113 , but 
we neglect these in our treatment. We do not have to 
concern ourselves with such contributions for A±, even 
in principle, since we can take the compression inside the 
lattice to be uniform in all directions. 

We also follow Fuchs in using the Ewald method to 
compute the resulting lattice sums. (See, e.g., |47| for a 
modern exposition of the Ewald method. Additionally, 
we make certain modifications to the method so that it is 
compatible with dimensional continuation — these are dis- 
cussed in Sec. lV Al ) In the Ewald method, one computes 
the sum of a slowly convergent or divergent series by ex- 
pressing it as the sum of two rapidly convergent series, 
introducing an Ewald screening function and using the 
Poisson summation formula. (Actually, the Poisson sum- 
mation formula is not strictly applicable to the functions 
being summed in many applications, which is what al- 
lows the Ewald method to regularize divergent series, by 
much the same method as zeta function regularization — 
see the Appendix for further discussion.) 

The Fuchs expressions are obtained by computing the 
energy of the perturbed lattice. Explicitly, we write the 
electrostatic energy per unit cell of the lattice (i.e., the 
energy required to remove a single blob) as 



W 



mm 

Q 



(14) 



[Compare the expression for the potential in Eq. (7) of 
Johnson and Ranganathan 14 71 and Fuchs's expression 
for the energy in Eq. (10) of [4^.] Here Q is the charge 
of a blob; A is the lattice under consideration, with dual 
lattice A*, and physical Voronoi cell volume We shall 
also use Conway and Sloane's [IH notation of VdetA 
for the Voronoi cell volume without the physical scaling. 
Primes on the sums denote the omission of the zero vec- 
tor; E is the Ewald screening function, with complement 
Ec{x) := 1 — E{x); <j) is the potential due to one of the 
blobs [including the effects of (physical) charge screen- 
ing] ; and the circumflex denotes the Fourier transform of 
(here three-dimensional) radial functions. (See Sec. IV Al 
for our Fourier transform conventions.) 



We now give the lattice sum expressions for the elastic 
constants. The most straightforward are those for cn and 
C44. which are of the form {(PW/ de'^)\^_^ /^l, for an ap- 
propriate deformation parametrized by e. (Note that we 
need cn in order to compute A±^ and the dimensionally 
continued Ajat, even though that elastic constant does 
not appear in /^off by itself.) The resulting expressions 
for elastic constants take the form (c e {011,044}) 



2^ 



E' 



{4>E)\\\x\\)^ + {4>Er(\\x\\)^^''^^ 



mrm) 



o^pW^ _ 2 dnd\\p\\ 

de^ de de 



de 

I — 



-y' 

d\\p\? ' 
de 



ueMp\\) fdn 



de 



(15) 



where the coordinate and 57 derivatives are given by 
Eqs. (Uni and ([TT]). [Recall that is unchanged by the C44 
perturbation.] The superscript bullets (•) denote deriva- 
tives taken with respect to — cf. the discussion below 

Eq. (nr 



requisite derivatives are thus 



d\\x\\ 



X1X2 



t=0 



de^ 



e=0 



To obtain the derivatives needed for the lattice sum ex- 
pression for C44, we note (following Fuchs ^3|) that C44 is 
the only elastic constant that contributes to the change 
of energy of the lattice if one considers a simple shear 
deformation. If we take this simple shear to be in the 
xi and X2 directions [now denoting our Cartesian coordi- 
nates by (cci, X2, cca), instead of the previous (x, y, z), for 
notational convenience], it corresponds to the deforma- 
tions xi — > x'l -|- ex2-, X2 ^ X2+ exi for the direct lattice, 
and pi {pi - ep2)/{l ~ e^), p2 -> (p2 - epi)/(l - e^) 
for the dual lattice (^3 and p^ are both held fixed) . The 



d\\p\r- 



de 



2piP2, 



de^ 



2pl 



(16a) 
(16b) 



[Nota bene (N.B.): We have given the derivatives of 
since the dimensionally continued Fourier transform de- 
pends naturally on this quantity — cf. Eq. p9)) . Addition- 
ally, all of the right-hand sides of these expressions are 
evaluated at e = (both here and in all similar situations 
in the next equation).] 
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For cii, we use the deformation xi — >■ (1 + e)xi, pi — >■ 
Pi/{1 + e) (with all other coordinates held fixed). This 
perturbation changes the cell volume, so we must include 
derivatives of 51, giving 



aiifii 



de 

d\\p\\' 



de 

on 
97 



E=0 



e=0 



9e2 



e=0 



e=0 



0. 



(17a) 
(17b) 

(17c) 



For Aiat, Fuchs uses the elongational shear in the xi 
and X2 directions, for which the change to the lattice's 
energy is given by 2Aiat alone (and has a lattice sum ex- 
pression similar to that for C44). However, as we shall see 
in Sec. |Vl the resulting expression is not well suited for 
dimensional continuation, since it contains fourth powers 
of more than one coordinate [see Eq. (12) in Fuchs [i^ ]. 
We can obtain an expression that is well-suited to di- 
mensional continuation if we note that Aiat = cn — C12 
(we have cn = C22 by cubic symmetry); Wc compute 
cii as above, but obtain C12 = 5*1122 + ^es by first us- 
ing a two-component perturbation to obtain S'1122, and 
then computing Pes o lo, Fuchs. The two-component per- 
turbation we use to obtain 6*1122 is a:i — > (1 + ei)a;i, 
X2 ^ (1 + e2).T2, Pi Pi/(1 + ei), P2 P2/(l + £2) 
(with 2:3 and p^ held fixed) , yielding the derivatives 



d\\x\\ 



dej 

d\\p\\ 



ei,2=0 
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de, 
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deide2 



,,2 ,,2 



, (18a) 



£1,2=0 



ei,2=0 



deide2 



n,2=0 11^11 = 

= 0, (18b) 



deide2 



ei.2=0 



(18c) 



where j S {1, 2} and ei_2 = ^ ei = £2 
obtains, from 5ii22 = {d'^W/d^id^2)\^^^ 
expression for W in Eq. 



= 0. One then 
P /ri and the 



5*1122 — 



E 



' xlxl 



mnwxw 



{mm\ 



-y' 
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pGA* 
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£1,2=0 



(19) 



We calculate P^^ using the same perturbation we used 
to obtain cn [with the derivatives given in Eq. (jl7p ]. 
Specifically, Pes = - {dW/de)\^^^ /fl [cf. Eq. recall- 



ing that Ski = -Pes4/]- This gives 



Pn. 



n 

peA* 
We thus have 



-y' 



(0i?,)(|b1|) + 2pf(0i?,)-(|bl) 



A 



lat 



cii - 5*1 



122 



Po. 



e=0 

(20) 



(21) 



where Cn is given by Eqs. and pT)) . and 6*1122 and 
Pes are given by the above expressions. 

Now, to obtain A\_ , we need to supplement the expres- 
sion for cii with the contributions due to changing the 
cell energy (since we are changing its radius) , along with 
the contributions due to changing the blobs' charge. The 
contribution due to changing the cell energy is 



A 



2 P, 



cell 



_L,coll 



(22) 



where Pccu/il is given in Eq. This can be deduced 
from the scalings of the Coulomb and surface contribu- 
tions [see Eq. (2) in Pethick and Potekhin [4l| ], noting 
that the change in the cell radius with this perturbation 
is given by 



dr 
'de 



r 



(23) 



which comes from noting that O = CcouT"'', where Cceii is 
a constant (since x is fixed, as we are keeping the overall 
density fixed) and (i9il/(?e)|c=o = ^ for the cn pertur- 
bation. For the derivatives of the blobs' charge, wc note 
that the (three-dimensional) charge density is fixed, so 
that the derivatives of the blobs' charge can be obtained 
from those of the cell volume [in Eq. (jl7cp ] by replacing 
the cell volume with the blobs' charge. Explicitly, we 
have 



dQ 
de 



de^ 



0. 



(24) 



We now show how to put together all these contribu- 
tions to obtain the sum that gives A±: We have 



Ai_ = cii + ^j_,ccii + Aj_^Q, 



(25) 



where cn is given by Eqs. p5|) and pT|) . A±^ccn is given 
in Eq. (|^, and 



xeA 



(,/)P)(||f||) + 2^(0P)'(||f|| 



i^' [{^Mp\\) + ^pii?E,n\\p\\) 



pGA* 



n J 



(26) 
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As one would expect from Earnshaw's theorem, cn is 
always negative, so one relies on the contributions from 
^±.ccii and A±^Q to make A± positive so that the lattice 
is stable to shears. (See Sec. IVII for further discussion.) 



V. DIMENSIONAL CONTINUATION OF 
LATTICE SUMS 

A. The dimensionally continued Ewald method 

In order to compute these lattice sums numerically, 
we employ a generalization of the standard Ewald [46[ 
method used by Fuchs [l^ . Showing the standard integer 
dimension version first, we have, summing a function / 
over a lattice A with dual lattice A*, 



1 



SeA 



VdcTA^^^. 



Y.ifE,)ip). (27) 



Here Ec{x) := 1 — E{x) is the complement of the Ewald 
screening function E, g{p) := g{x)e~^'^^^'^d^x de- 
notes the standard Fourier transform, and we choose E 
so that both of the sums on the right-hand side converge 
quickly. [Recall that the classical Poisson summation for- 
mula says that 'EseA fi^) = (det A)-!/^ J^peA- fip)-] 

The classic choice for E for Coulombic potentials (dat- 
ing back to Ewald) is the complementary error function. 
However, this turns out to be insufficiently flexible to pro- 
vide good convergence for the sums we consider (partic- 
ularly for small d). Following Nijboer and de Wette [s^ 
and Fortuin [stI ]. we introduce the incomplete gamma 
function, r(-, •), and use the screening function 



E{x)^r{N/2,aW)/mm. 



(28) 



This reduces to Ewald's complementary error function 
for iV = 1. The extra freedom contained in N allows us 
to tune E to provide fast convergence for the sums we 
encounter. We used iV = 10 and a = 1.2/ a in the compu- 
tations reported in Sec. lVIl [Here a is the lattice spacing, 
given in Eq. (j^Hl)-] We could have doubtless obtained 
faster convergence if we had allowed these parameters 
to vary with d (and possibly also A), particularly for d 
close to 1. However, we found that these values to give 
reasonably good performance, and thus did not perform 
much experimentation beyond checking that our results 
are insensitive to small variations in the Ewald screening 
parameters. 

In order to dimensionally continue our lattice sums, we 
need to dimensionally continue the Poisson summation 
formula. We shall give an overview of the calculational 
aspects here — see [SOj for more details of the derivations, 
and a proof of the formula for nicely behaved functions. 
We first introduce the dimensionally continued Fourier 
transform for spherically symmetric functions, given in 
Theorem 3.3 of Chap. IV of Stein and Weiss 



9{p) := 



2tt 



P 



,d/2-l 



gir)Jd/2-ii^7rpry^^dr. (29) 



Here is the Bessel function of the first kind. (See 
Sec. 2.2 of [s^l for more details, including an alternative 
expression in terms of a hypergeometric function.) We 
also note that we can compute the dimensionally contin- 
ued Fourier transform of the potential from its defining 
partial differential equation (l8|), giving 



47r 



(30) 



[This agrees with the result of the more involved calcula- 
tion one could perform using the dimensionally continued 
Fourier integral given in Eq. ([^ .] We use Eq. (PH)) along 
with the dimensionally continued Fourier integral (|29p 
to compute (j)Ec efficiently, writing it as ^ — (j)E, where 
the integral giving the second term converges reasonably 
rapidly. 

We then introduce the theta series of a lattice. As is 
discussed in more detail in Sec. 2.3 of Chap. 2 of Con- 
way and Sloanc j5l[, the theta series of a lattice is the 
generating function of the number of lattice points on a 
sphere of a given squared radius, so the theta series of a 
lattice A is defined by 



llfell 



(31) 



feeA 



[N.B.: There are several different notational conventions 
for theta series and theta functions. We have chosen to 
write all our theta series and functions as functions of the 
nome, q, unless we note otherwise. These exceptions will 
only occur in discussions of the Jacobi formula for the 
theta series of the dual lattice, where it is convenient to 
treat the theta series as a function of a complex variable 
z, with q ~ e*'^^. We shall denote this by an overbar — 
e.g., Qa{z) ■= 9A(e*'^^). Conway and Sloane treat all 
their theta functions as functions of z, even when they 
write their expansions in terms of the nome.] Thus, if 
we define the power series coefficients of the theta scries 
using 



(32) 



we can write the sum of a spherically symmetric function 
F over A as 



(33) 



k£A 



1=0 



The dimensionally continued Poisson summation for- 
mula then has the form 

oo oo 

E AiFi^i) = -== E A^HVBf)- (34) 
J^o vdctA^ 

Here, A^* and _Bj* are defined analogously to their coun- 
terparts for A. The theta scries for A* can be calcu- 
lated from 8a using the Jacobi formula ([M)) . which is 
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already in dimensionally continued form. [N.B.: The di- 
mensionally continued Poisson summation formula pre- 
sented in [501 omits the factor of (dctA)^^/^, since, as 
discussed there, this factor cancels against a similar one 
present in the Jacobi transformation formula. We in- 
clude the factor here, since we obtain a simpler result 
for the theta series of the dual lattice by using the stan- 
dard Jacobi transformation formula.] Of course, we need 
to sum more than just spherically symmetric functions 
[see Eqs. (fT5 |) -(f20 | ]. but, as we shall see shortly, this di- 
mensionally continued Poisson summation formula will 
be sufficient for our needs. 



B. Dimensional continuation of the lattice 

Turning now to the problem of dimensionally contin- 
uing the lattices themselves, recall that the integer di- 
mension lattices are [up to an overall scaling, which we 
determine in Eq. (|49|) ] a bcc lattice for d = 3, a hexago- 
nal lattice for = 2. and Z, the (one-dimensional) lattice 
of integers, for d = 1. In some sense, the natural way to 
dimensionally continue these lattices would be to dimen- 
sionally continue the root lattice family A'^, which gives 
those lattices for d G {1,2,3} and, more generally, gives 
the best lattice covering of R" for rt < 5 (as discussed 
in Conway and Sloane [5l|). However, unlike most other 
families of lattices, does not have a theta scries that 
is written in a nicely dimensionally continued form. Its 
theta series is written in terms of a sum whose number 
of terms depends upon dimension — see, e.g., Eq. (56) in 
Chap. 4 of Conway and Sloane [5l| for the theta series 
for Ad, from which the theta series for the dual lattice 
can be obtained using Jacobi's formula [our Eq. ([39])]. 
The sum is over dth roots of unity, so one could contem- 
plate writing it as an integral, using Cauchy's theorem, 
and proceeding that way. However, even disregarding the 
complications this would involve, it is not clear how to 
compute the sums involving xf that we need (for, e.g., 
cii) in this framework, or, alternatively, how to imple- 
ment the requisite distortions to the lattice at the level 
of its theta series. 

We thus proceed by treating the dimensionally contin- 
ued lattice as a union of hyperlattices (i.e., lattices of 
one dimension fewer than the overall lattice)^ whose sep- 
aration is given by a freely specifiable function /^^* that 
interpolates between the integer dimension separations. 
One then finds that the sums over the hyperlattices di- 
mensionally continue in a natural way, and that the final 
result for the shear modulus is rather insensitive to the 
choice of Z'*^* (provided that it satisfies some reasonable 
properties, discussed below). 

Explicitly, the (scaled) integer dimension lattices can 



be written as [/'=*' (d)Zovcn] x Z'*"! + [f''*'{d)Zodd\ x 
[Z"^^^ -I- (1/2)''"^], where Zoven and Zodd denote the even 
and odd integers respectively, denotes the (d — 1)- 

dimensional lattice of integers, and Z"*"^ -I- (1/2)'^"^ de- 
notes the same shifted by the [{d— l)-dimensional] vector 
all of whose components arc 1/2. This is illustrated in 
Fig. m where we take xi to be in the direction orthogonal 
to the hyperlattices. 





• 
■ 


1 Xl 




Z^-i + (l/2)''-i\^ I 


















z^-'-^ I 



FIG. 4: A schematic of the decomposition of the full lattice 
into hyperlattices. 

Here /'^'(d) is freely specifiable, except that it must 
satisfy {1, 2, 3} ^ {1, a/3/2, 1/2}. Additionally, it makes 
sense to choose it to be smooth, nonincreasing, and con- 
cave. A possibility that satisfies all these criteria is 



nilid) :=cos(^-[d-l]). 



(35) 



which we will use in all our results, unless otherwise 
noted. We will show that the results for other possi- 
bilities all agree well with /^q*. Consider the envelope of 
all possible functions satisfying the requirements. The 
bounds on this envelope are 



_UV3/2-l){d-l) 
\[(l-V3)(d~2)-f 



f 1, 
731/2, 



[(1- V3)(d-2)H 
(x/3/2-l)(d-l) 



V3]/2, 
+ 1, 



1 < d < 2, 

2 < d < 3, 



(36a) 

l<d<d', 
d' <d<2, 
2 < d < 3, 

(36b) 



Note that many of these hyperlattices arc actually shifted lat- 
tices, mathematically speaking. 



where d' := 3 - 1/(^3 - 1) ~ 1.63. See Fig. [5] for an 
illustration. While it would be natural to choose Z''^' so 
that the dimensionally continued lattice had effective cu- 
bic symmetry (discussed in Sec. IV C[) . it is not possible 
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to do so. If we consider, for instance, d = 5/2, then the 
Z'^~^ hyper lattices have theta series l + 3q + 0{q^), while 
the theta series for the Z'^^^ + (1/2)''"^ hyperlattices is 
2'^/^g'^/^[H-0(g^)]. The irrational prcfactor in the second 
theta series is a clear obstruction to obtaining effective 
cubic symmetry: There is no way of obtaining an irra- 
tional contribution from the first theta scries, since all its 
coefficients are rational. 




FIG. 5: Plots of /^S [ given in Eq. (|35|l ]. along with the point- 
wise sup and inf over all such possibilities [/iup and //,jf , given 
in Eqs. p6l)]. 

Now, we dimensionally continue the sums over the hy- 
perlattices using the hyperlattice's theta series, which are 
naturally dimensionally continued. Specifically, the theta 
series of Z'^ and + (1/2)'' are i?f and where 

Mq) E'?^'^'^'^'' ^= E^'". (37a) 

fcez fcez 

fcez 

See, e.g.. Sec. 5 of Chap. 4 of Conway and Sloane [sij . 
but recall the differences in their theta function notation, 
discussed below Eq. (PT|) . We also introduce here since 
it appears in the theta series of the dual lattice. 
The theta series of the full lattice is thus 

= ^3([2/'^'(rf)]'^) [^3(^)]''' 

One can check that this reduces to the appropriate ex- 
pressions for d G {1,2,3}; the theta series for the 3- 
dimcnsional bcc lattice and 2-dimensional hexagonal lat- 
tice are given in Eqs . (96) and (60) of Chap. 4 of Con- 
way and Sloane [51[. For d = 1, we obtain the theta 



series for Z in a nonstandard form — one can convert it 
to the standard one given above by using the identity 
z?2(4z) + i?3(4z) = Mz) from Eq. (22) in Chap. 4 of 
Conway and Sloane. 

We use Jacobi's formula (which is already dimension- 
ally continued) to obtain the theta series of the dual lat- 
tice, viz., [e.g., Eq. (4) in [50j . or Eq. (19) in Chap. 4 of 
Conway and Sloane [5l[ ] 

Qa,{z) = VdctA{i/zf/^eA{-l/z). (39) 

When applying Jacobi's formula, we take the volume of 
the lattice's Voronoi cell to be dimensionally continued in 
the obvious way, viz., Vdet A = /''''(d). [N.B.: This ex- 
pression neglects the overall scaling of the lattice, which 
we fix in Eq. (|49|) .] We then obtain, upon use of the theta 
function identities in Eq. (21) of Chap. 4 in Conway and 
Sloane, 

e^'iz) ^ {Mz/[2f^\d)r)[Mz)Y-' 

([^3(,)]-V[^4(.)]'^-^) 
^^[(fe+l/2)//-(rf)]^ ([^^3(^)]'"'- [^4(Z)]'"')}. 

(40) 

The second expression tells us that we can treat the dual 
lattice as a union of hyperlattices separated by a distance 
of [2/'^'((i)]^^, where the even hyperlattices are Dd-i 
and the odd ones arc Dd-i + (O'^^^l) (cf. the schematic 
of the direct lattice shown in Fig. |4|). Here Dd is the d- 
dimensional root lattice discussed in Sec. 7.1 of Chap. 4 
of Conway and Sloane [Slj, and Dd + (0'^~^1) denotes 
the same lattice shifted by a unit in one coordinate di- 
rection. These lattices' theta series are [Eqs. (87) and 
(89) of Chap. 4 of Conway and Sloane ^ ] (i^^ ± i?|)/2, 
with the upper [resp. lower] sign corresponding to Dd 
[resp. Dd + (0'^~'^1)]. This interpretation is convenient, 
as it allows us to compute the sums over the dual lattice 
using the same technology as for the direct lattice. 

C. Lattice sums 

From the expressions in Eqs. p5|) "p0 | . we see that we 
need to dimensionally continue sums involving a spheri- 
cally symmetric function times either x" (n € {0,2,4}), 
or x\x\. In order to do this with our method, we note 
that we can express x" in terms of f^^^{d) and the index 
of the hyperlattice (fc, in our discussion below), due to 
our method of dimensionally continuing the lattice. For 
X2, we note that the hyperlattices have effective cubic 
symmetry, so we have 
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where Hd-if^Sr denotes the intersection of a hyperlattice 
Jid-i (with dimension d — I) with a (d — 2)-sphere of 
radius r centered at the origin, and n-fi (r) is the number 
of points in this intersection (obtained from Hd-i's theta 
series). Unfortunately, this method is not applicable to 
X2 (which appears in the Fuchs expression for yliat along 
with xf), which is why wc compute Aiat using the more 
involved method discussed above. However, note that we 
can compute Pes (and thus Ajat) and C44 two ways — the 
lattice sums we use do not contain fourth powers of any 
of the coordinates, so we can also compute them with the 
index substitution 1 o- 2. We find that the two choices 
differ by < 10%. Since these would be exactly equal if 
effective cubic symmetry held, this gives an indication 
that our assumption of such symmetry is valid to about 
this level. 

Wc now give the explicit expressions for the computa- 
tion of the requisite lattice sums in a dimensionally con- 
tinued manner. If we consider summing G{xi, 2:2)^(11x11) 
over the d-dimensional lattice A^, where G{xi,X2) is ei- 
ther one of cc" (n € {0,2,4}), cc^, or x^x"^, then we have 
to evaluate the following double sum: 



[recall that the hyperlattices' theta series are d^-^ and 
i?3~^; the theta functions are defined in Eq. p7ap ]. Sim- 
ilarly, for the dual lattice, we have 

00 00 

'=0 (46) 
where 

Cpf = [k/f^\d)r, = [(fe + 1/2)/ f^\d)]\ 

(47a) 

=2//(d-l), V„. = i2l + l)/{d-l), (47b) 



00 00 

^ G(xi,.T2)F(||x||) =^^{.,7V/^UGF(7^[31) 

fc=0 1=0 



(42) 



where 



= [2kf-\d)r, = [{2k + l)/^*(d)]", (43a) 

A,2^l/id-l), B,2=l/id^l) + l/A, (43b) 



7^[3l := ^[2kf^\d)]^ +1, (43c) 
Tel'l := y^[(2fc + l)/i'^t(d)]2 + / + (d-l)/4, (43d) 
and one obtains the coefficients for the 

■^1*^2 ^^^^ 

by mul- 
tiplication (e.g., Ax2y2 = Ax^Ayi). [The superscripts 
[2] and [3] come from the names of the theta functions 
that give the pertinent hyperlattices' theta series — see 
Eq. (|45p .] The Z-sum comes from summing over an in- 
dividual hyperlattice, and the fc-sum then sums over all 
hyperlattices. We have introduced 



1, fc = 0, 



otherwise, 



(44) 



so we can take k to run only over the positive integers, 
while the hyperlattice index runs over all integers. The 
structure of the hyperlattices is accounted for by the 
theta scries coefficients N; , defined by 



[3] J 



/=0 



Ml) 
0I/4 



(45) 



k=0 1=0 

+ 2Nr 



TZ+ :=^[k/f-^{dW + 2l, (47c) 
7^- y^[(fc+l/2)//i-t(d)]2 + 2/ + l, (47d) 

(with the same multiplication for the P1P2 case as for the 
direct lattice), and the theta series coeSicients are given 
by 



(48) 



[recall that the dual hyperlattices' theta series are iJg ^ ± 

We can now calculate the elastic constants by combin- 
ing together these results with our previously derived ex- 
pressions. For C44, these are given by Eqs. ([T5|) and (fT6|) . 
For ^lat and A±, one uses the sums given in Eqs. ([2T|) 
and (pS)) . respectively. We also need to determine the 
scaling of the lattice for a given energy density. This is 
given by the cell radius, R, computed in Sec. Ill Cl Equat- 
ing the volume of a cell with this radius to the volume of 
the lattice's Voronoi cell [given by /'^'(d)a''] determines 
the lattice's spacing, a. Explicitly, we have 



rl/2 



[P^^{d)T{d/2 + l)]^/d 



(49) 



using the dimensionally continued expression for the vol- 
ume of a unit d-baU, viz., n'^/'^ /T{d/2 + 1). Additionally, 
we compute Q, the charge of the blob, using the charge 
density of the rare phase and the (d-dimcnsional) volume 
of the blob. 



VI. RESULTS AND DISCUSSION 

Here we present the shear moduli for the various EOS 
parameters we consider (given in Table |l|. 
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• Dimensionally continued 
a No dimensional continuation 



Max. pressure in stable OV stars 




maximum mass star in Fig. [T) 



2 3 4 

Pressure (10''^ erg cm 



FIG. 6: The effective shear modulus and lattice dimensional- 
ity versus pressure for the Hyl EOS both with and without 
dimensional continuation. We have also indicated the maxi- 
mum pressure one obtains in a stable OV star using this EOS. 




FIG. 7: The effective shear modulus and lattice dimensional- 
ity versus OV (Schwarzschild coordinate) radius for the max- 
imum mass stable star (of total radius 12.5 km) with the Hyl 
EOS both with and without dimensional continuation. 



First, to give an indication of the effects of dimensional 
continuation, in Fig. [5] we plot the effective shear modu- 
lus /ieff [see Eq. (|13p] and lattice dimensionality d versus 
pressure with and without dimensional continuation, for 
the Hyl EOS with a surface tension of ct = 80 MeV fm^^ 
[and the /^q' lattice interpolation function from Eq. ([55]) ]. 
Wc have shown all the pasta phases, even though only 
the first few appear in stable stars, as is indicated in 
the figure. The jump in the effective shear modulus oc- 
curs at the halfway point — i.e., equal amounts of quark 
and hadronic matter — and is due to the switch on the 
screening length A at that point (see the discussion in 
Sec. IIID[) . Also note that the lattice becomes unstable 
for a small range of d slightly less than 3 — see the discus- 
sion below. Wc additionally show fics and d versus the 
(Schwarzschild coordinate) OV stellar radius Tgtar for the 
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FIG. 8; The (dimensionally continued) Aiat versus pressure 
for the Hyl EOS and a variety of /'"^s. See the text for a 
discussion of the small unphysical dip below zero. 
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FIG. 9: Ratios of the effective shear modulus computed us- 
ing /sup and /i'^f to that computed using flU [see Eqs. (PS]) 



and (|36p ] plotted versus pressure for the Hyl EOS and a = 
80 MeV fm-^ Thus the error due to our lack of knowledge 
of this function is at most about 5% in small regions of the 
star. 

Even if the effective shear modulus is positive, the lat- 
tice can still be unstable to shear strains if ^lat is neg- 
ative. This is the case for small regions of the lattice 
whose effective shear modulus is shown in Fig. [6l specif- 
ically, where d is slightly less than 3. However, one only 
obtains such an instability when one uses /^q*, //^j, or 
something similar for the lattice interpolation function. 
If one uses /gj^p, then ^lat remains positive. This is il- 
lustrated in Fig. El for the same Hyl a = 80 MeV fm"^ 
case, using the lattice interpolation function fl^l as well 
as the sup and inf over all interpolation functions [given 
in Eqs. (|36p ]. We do not use fl^^ itself in our shear mod- 
ulus calculations because it docs not satisfy the smooth- 
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ness and convexity criteria, though Fig. [8] indicates that 
an /'^' that was somewhat greater than fl^l for d near 3 
would not lead to the instabilities. Regardless, the varia- 
tions in shear modulus are at most 5% between sup and 
inf in small regions, as illustrated in Fig.|9l so we tolerate 
this small inconsistency. 

The sign issue does not arise with which is made 
positive by the contributions from changing the cell size 
[cf. Eq. (f25|) ]. While cn is always negative, A±^cc\\ is 
enough to make A± positive except for d ~ 2.99, where 
A±^Q is also needed in certain cases. There A±^q is pos- 
itive, although it is negative for d <2. 
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FIG. 10: The (dimensionally continued) effective shear mod- 
ulus versus pressure for the Hyl EOS and a variety of surface 
tensions, a. The discontinuity, which does not occur in a real 
star, is largely caused by our simplified treatment of charge 
screening and serves as an estimate of the error in that treat- 
ment. For high or moderate surface tensions, that error is less 
than a factor 2. For low surface tension the error can be an 
order of magnitude, but these lattices are nearly unstable in 
any case, as shown below. 
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FIG. 11: The (dimensionally continued) effective shear mod- 
ulus versus pressure for the Hyl EOS and two choices of a, 
showing the negative values for a = 10 MeV fm""^. 



The largest effect of the equation of state parameters 
on observable quantities involving the shear modulus is 
the extent of the mixed phase in stable stars. See Weis- 
senborn et al. [1^ for a survey of the dependence of the 
mixed phase's extent on the bag constant and QCD cou- 
pling constant for two hadronic EOS parameter sets. 

The largest effect on the shear modulus itself, however, 
comes from the surface tension, a. This is illustrated for 
the Hyl EOS in Fig. [TUl (The magnitude of the reduc- 
tion of the shear modulus with decreasing surface tension 
also holds for the other EOSs we consider.) Note that 
the shear modulus is greater for smaller surface tensions 
for d close to 3. Additionally, a has a direct, and quite 
substantial, effect on the lattice's stability: If a is too 
small, then the lattice will become unstable to shears (in- 
dicated by A± becoming negative) as the dimensionality 
decreases, as illustrated for the Hyl EOS and two small 
values of a in Fig. [TT] A surface tension that is much 
higher than the range we consider makes any lattice too 
energetically expensive to form. 

To illustrate the relatively small effect of the other EOS 
parameters on the effective shear modulus for a fixed sur- 
face tension, we plot ^eS versus the quark volume fraction 
X for a representative sample of the EOSs from Table |T] 
with cr = 80 MeV fm"^ in Fig. [H The different fla- 
vors of Hyl EOS— i.e., Hyl, Hyl^, Hylcr, Hyl^cr— all 
have quite similar shear moduli for a fixed surface ten- 
sion. The inclusion of the surface tension contribution to 
the pressure balance causes the largest difference of any 
of these other EOS parameters, but even this difference 
is considerably smaller than the differences between the 
EOSs shown in the figure, and is negligible where the 
shear modulus is largest. We have thus not shown the 
traces for those EOSs, to avoid a cluttered figure. 




0.4 0.6 
X 



FIG. 12: The (dimensionally continued) effective shear mod- 
ulus versus x for a representative selection of the EOSs from 
Table Hand ct = 80 MeV fm~^ We have left off a few points 
at either extreme of x with low shear moduli to better show 
the differences between the different predictions where the 
shear moduli are the highest and most astrophysically rele- 
vant. There the uncertainty in shear modulus due to factors 
other than surface tension is at most about a factor of 2. 
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It is interesting to compare these shear modulus val- 
ues to those for the lattice of nuclei in the neutron star 
crust as well as those for crystalline color superconduct- 
ing (CSC) quark matter: The crustal shear modulus 
ranges from ~ 6 x 10^^ to ~ 2 x 10'^° erg cm~'^ (see, e.g., 
Fig. 2 in (111), while the shear modulus of CSC quark 
matter computed by MannarcUi, Rajagopal, and Sharma 
(MRS) could be as large as ^ 4 x 10 erg cm ^; 
the lower bound they give is ^ 8 x 10^^ erg cm~^. The 
hadron-quark mixed phase shear moduli we have com- 
puted range from ^ 10'^° erg cm~^ for small blobs to 
~ 8 X 10'^'^ erg cm~'^ for hadronic slabs (with a surface 
tension of ct = 80 MeV fm~^). Thus the mixed-phase 
shear modulus is at least comparable to the largest shear 
modulus in the crust and can be three orders of magni- 
tude higher than it, and the largest mixed-phase shear 
modulus is an order of magnitude less than the largest 
CSC quark matter shear modulus, roughly the geometric 
mean of the range given by MRS. (Note, however, that 
the maximum mixed-phase shear modulus for a given set 
of EOS parameters is not necessarily realized in a stable 
star.) 



VII. CONCLUSIONS 

We have made a more careful calculation of the shear 
modulus of the hadron-quark mixed phase than has pre- 
viously been attempted [H, H^, . In particular, we 
have computed all the lattice's (anisotropic) shear elastic 
constants, before averaging to obtain an isotropic effec- 
tive shear modulus for a polycrystal. We have also dealt 
with the lattice's changing dimension in both the electro- 
static potential and the geometrical effects on the elastic 
constants. Perhaps most importantly, we have included 
the contributions to the elastic constants from changing 
the size of the blobs for d < 3. These act to stabilize 
the lattice for lower dimensions (leading to significant 
contributions to the shear modulus from these portions 
of the lattice), though only for sufficiently large surface 
tensions. We have found that for our choices of parame- 
ters, most of the lower-dimensional portions of the lattice 
are unstable to shear perturbations if the surface tension 
is less than 10-20 MeVfm-^. (As discussed in [Ullll, 
the mixed phase is not favored if the surface tension is 
too large. However, as we mention in Sec. Ill CI and will 
explore in depth in (36j . these are local energy consid- 
erations, while the mixed phase may still be favored by 
global energy considerations.) 

These calculations depend upon a wide variety of 
poorly constrained parameters. While we found that the 
shear modulus of the mixed phase itself depends most 
sensitively on the surface tension, astrophysical effects 
depend on the amount of mixed phase present in a given 
star, which is primarily determined by the standard hy- 
brid EOS parameters. We will see this in more detail 
when we use the shear moduli calculated here to com- 
pute the maximum clastic quadrupolar deformations of 



these hybrid star models in a companion paper |36| . 

The obvious place where the shear modulus calculation 
could be improved significantly is the treatment of charge 
screening: We know from the calculations of Endo et 
al. [40I [ooj that including (nonlinear) charge screening in 
the computation of the cell energy leads to significant dif- 
ferences in the lattice properties at a given density. Here 
one would like to at least use the Thomas-Fermi approx- 
imation, if one did not perform a nonlinear calculation 
as in [io, [13] ■ The shear modulus will be affected both 
by the change in cell size and spacing, as well as through 
the cell energy's direct contribution to one of the elastic 
constants for lower dimensions. Since this contribution is 
necessary to stabilize the lattice for those dimensions, it 
would be interesting to see whether our discovery that the 
lattice is only stable for sufficiently large surface tensions 
still holds with charge screening. For instance, there will 
be, in effect, Fermi contributions to the cell energy in a 
proper treatment. 

One might also want to investigate the effects of using 
different descriptions of the hadronic and quark matter. 
For instance, for the quark matter, one could use the 
higher-order perturbative calculations of Kurkela, Ro- 
matschke, and Vuorincn [6^ or the Nambu-Jona-Lasinio 
treatment (e.g., [.22.,„91j), while Dirac-Brueckner-Hartree- 
Fock (e.g., [91|) would be a possibility for the hadronic 
matter. One might also investigate the potential ef- 
fects of color superconductivity (see Q for a review): 
The inclusion of CSC quark matter would, of course, 
increase the shear modulus (see Mannarelli, Rajagopal, 
and Sharma for a calculation of the shear modulus 
of bulk CSC quark matter). However, if one does not 
have a crystalline phase, color superconductivity would 
primarily affect these shear modulus calculations through 
the EOS, where Alford et al. [65| find that including color 
superconductivity reduces the transition density to quark 
matter, but does not appreciably change the maximum 
mass. While one would also expect color superconductiv- 
ity to affect the energy of the quark blobs, and thus the 
contributions to the shear modulus from changing the 
blob size, we have not even included quark interaction 
contributions to the blob energy in the present calcu- 
lation, only considering the electrostatic contributions. 
Other possibilities for extending the calculation overall 
would be including more exotica (hyperons, for instance, 
as in, e.g., [Ill, [22|, [5^ ) , and, in particular, magnetic fields 
[cf. the discussion in Baiko [i^ and below our Eq. ([T^ ]. 

One could also apply our methods of calculating the 
shear modulus of the pasta phases to the nuclear pasta 
appearing in the crust. Such a calculation would be par- 
ticularly interesting given recent studies of the observa- 
tional effects of the crustal pasta which used very rough 
models for its shear modulus [H, H^l . Our methods could 
also be adapted to calculate the shear moduli of me- 
son condensates more carefully than the existing order of 
magnitude estimates (Ref. [2a] and references therein). 

The shear moduli here can be several hundred times 



larger than those first estimated [12|] . Naively, one might 
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expect the maximum quadrupole to go up by a similar 
factor, but this is comphcated by issues of where the 
various lattices occur in the star. We will present an 
exploration of that and other issues with the maximum 
quadrupole elsewhere [36| . 
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Appendix: Checks of lattice sums 

We have checked that our code can reproduce all the 
relevant results for elastic constants presented in the 
literature, and detail these checks here: For a three- 
dimensional unscreened bcc lattice, we have checked that 
we can reproduce the results of Fuchs [i^ (for C44 and 
^lat) and Ho [11] (for cn), and also that we are in agree- 
ment with the very recent and much more precise cal- 
culation of all three elastic constants (to 8 significant 
digits) by Baiko [i^. For a three-dimensional bcc lat- 
tice with screening, we have checked against the results 
from Horowitz and Hughto [95|. We have also checked 
that the surprisingly simple results we obtain in the 
unscreened two-dimensional hexagonal case agree with 
those obtained analytically using zeta function regular- 
ization. 

Numerical agreement is good. We reproduce the Fuchs 
results to better than 0.25% for ^lat and 0.015% for C44. 
Moreover, we agree with Baiko's results to all five signif- 
icant figures to which we have calculated them (using a 
screening length of lO^^a). Baiko claims agreement with 
Fuchs and does not comment on the discrepancy, which 
wc conjecture is due to modern computational technol- 
ogy allowing for the summation of more terms. We agree 
with Ho's result to the four significant figures to which 
he gives it. Horowitz and Hughto computed A\at (twice 
their 611) for a screening length of A = 0.863a. We agree 
to better than 1% for Ai^t and 0.011% for C44. The dis- 
crepancy for C44 could be due to rounding. The discrep- 
ancy for Aiat is compatible with Horowitz and Hughto's 
caveats about a 1% error in that part of their calculation 
due to using a brute force summation rather than the 
Ewald method.^ 

For the two-dimensional unscreened hexagonal lattice, 
we obtain what appear numerically to be simple fractions 



N.B.: Horowitz and Hughto's expression for 611 in terms of cn 
and C12 in their Eq. (13) is missing a factor of 1/2 , which can 
be seen to be present by starting from Eq. (7) in |83| and noting 
that C31 = C12 for a cubic lattice. 



for the elastic constants (when expressed in Fuchs-type 
units of Q^)- Explicitly, we have Aiat = 1/2, cn = —1/4, 
and C44 = 1/4. 

We can show that these clastic constants indeed take 
the expected, extremely simple forms by calculating them 
analytically using analytic continuation of an Epstein 
zeta function, as discussed in, e.g., [sl, 113] • (A cal- 
culation of these elastic constants for a 2-dimensional 
Coulomb lattice does not appear to exist in the litera- 
ture.) We first note that the 2-dimensional Coulomb po- 
tential is —2(3 log r. Thus the sums that give the elastic 
constants will all have the formal (and divergent) form 



a, 



(A.l) 



where a is some constant and Ahox is the hexagonal lat- 
tice. To see this explicitly for cn, we have (writing it in 
Fuchs-type units) 



hex _ 



- 2- 



where we have used the substitutions xf — )• ||x|p/2 and 
xf — > (3/8)||x||'* (valid when summing over a hexagonal 
lattice). (One can obtain these substitutions by direct 
calculation, for which it is convenient to write the lattice 
points as elements of C, so one can multiply a given point 
by sixth roots of unity to obtain all the other lattice 
points the same distance from the origin.) 

Of course, this sum diverges, but, as indicated in 
Sec IIVI we expect to obtain a regularized version from 
our Ewald sum procedure. Indeed, it seems reasonable 
(possibly even likely) that we should obtain the same re- 
sult using Ewald sums as one obtains using zeta function 
regularization, since both methods rely at their root on 
the Jacobi imaginary transformation of theta functions: 
See, for instance, the discussion in [50[, which obtains the 
Poisson summation formula used in the Ewald method 
from the Jacobi transformation, and the discussion in Ap- 
pendix B of [96j , which analytically continues the Epstein 
zeta function using the Jacobi imaginary transformation 
of theta functions (referred to by its alternate name of 
the modular property of the theta function). 

We thus begin by defining the Epstein zeta function 
Ca associated with an arbitrary lattice A. This is given 
[cf. the version for Z,'^ in Eq. (Bl) of [9^ and the more 
general version in Eq. (2.5) of [l^l ] by 



Ca(.) := E'' 



1 



2s 



(A.3) 



SeA 



for Re s > d/2 (where d is the dimension of A), and ana- 
lytically continued to s G C (except for isolated poles) 
by Riemann's integral-splitting procedure and the Ja- 
cobi imaginary transformation, as in [96l . l97j . We have 
Ca(0) = — 1 for any lattice [of positive dimension — cf. the 
third property given in Appendix B of [96| and Eq. (2.13) 
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97| ], so we obtain cn = —1/4, in Fuchs-type units, need to use the additional hexagonal lattice sum substi- 



as advertised above. We obtain the values for C44 and A tution 
given above by the same procedure, except that we also 



X2 ^ ll.fr ; 
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